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Abstract. The goal of this paper is to provide an analysis of the "toolkit" method used in the 
numerical approximation of the time-dependent Schrodinger equation. The "toolkit" method 
is based on prccomputation of elementary propagators and was seen to be very efficient in the 
optimal control framework. Our analysis shows that this method provides better results than 
the second order Strang operator splitting. In addition, we present two improvements of the 
method in the limit of low and large intensity control fields. 
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1. Introduction 

The control of the evolution of molecular systems at the quantum level has been a long 
standing goal ever since the beginning of the laser technology. After an initial slowed down of 
the investigations in this area due to unsuccessful experiments, the realization that the problem 
can be recast and attacked with the tools of (optimal) control theory [1] greatly contributed to the 
first positive experimental results [H [31 IH El [5] . Ever since, the desire to understand theoretically 
how the laser acts to control the molecule lead the investigators to resort to numerical simulations 
which require repeated resolution of the Time Dependent Schrodinger Equation of the type H]) ; 
additional motivation comes from related contexts (online identification algorithms, learning 
algorithms, quantum computing [7], etc.). 

The numerical method used to solve the time dependent Schrodinger equation must provide 
accurate results without prohibitive computational cost. The conservation of the LF' norm of 
the wave function ip(x,t) is also generally required for stability and as a mean of qualitative 
validation of the numerical solution. 

In this context, the second order Strang operator splitting is often considered [H |9l [T0| . 
However, this method suffers from two drawbacks. First, the numerical error is proportional to 
the norm of the control which implies poor accuracy when dealing with large laser fields e{t) and 
make necessary the use of small time steps. Secondly, it requires at each time step three matrix 
products. This difficulty is enhanced in some particular settings e.g., in optimal control, where 
the matrices involved in the control term must be assembled online. 

Recently introduced, the "toolkit method" [11] [12] solves this last problem by precomputing 
a set of elementary matrices, used in the numerical resolution. Each matrix is associated to (one 
or several) field values and enables to solve the evolution over one time step. This algorithm 
has been used in various frameworks and shows excellent results. It has also been coupled 
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successfully with optimal control and identification issues [13]. The dependence on the L°°- 
norni of the control, which is a restriction of the Strang method, is also improved by the "toolkit 
method" as it will be shown in our analysis. 

The goal of the paper is to provide a (first) numerical analysis of the "toolkit method". 
Our mathematical tools arc related to that in [8] (but for a different setting; see also [Ml [15] 
for connected results); the treatment here is different because of the quantization appearing in 
the values of the control s{t) which impacts both the mathematical analysis and the numerical 
efficiency of the method. The analysis enables us to propose two possible improvements. 

The paper is organized as follows: after having introduced the model and some notations in 
Section [2] the "toolkit" method is presented and analyzed in Section [3] An improvement of this 
method in the limit of small control fields is introduced in Section ID A second improvement, in 
the limit of large control fields is given in Section [5j Finally, Section [6] gathers some numerical 
results. 



2. Model and notations 

In this section, we present the Schrodinger Equation that will be considered in the paper and 
some useful notations. Note that the algorithms we consider in this paper actually apply for 
other types of Schrodinger equations. Indeed, they are rather based on the algebraic structure 
of the control problem (bilinear control) than on the regularity of the solution. In this way, the 
error estimates that we will obtain hold not only in L^, but also in H^. 

We consider the time dependent Schrodinger equation (TDSE) (7 G N): 

.^s ( idti^{x,t) = {Ho- n{x)eit))^{x,t), xeW 

^ ' \ il'ix,0) = ipoix), X e Wi. 

This equation governs the evolution of a quantum system, described by its wave function -0, 
that interacts with a laser pulse of amplitude e, the control variable. The factor fj, is the dipole 
moment operator of the system. The Hamiltonian of the system is Ho = + V where A^^ 
is the Laplacian operator over the space variables and V = V{x) the electrostatic potential in 
which the system evolves. We refer to [16] for more details about models involved in quantum 
control. Note that to obtain Eq. (jlj , one has considered the laser effect as a perturbative term, 
so that the control term e{t)fj,{x) is obtained through a first order approximation with respect 
to e{t). While often considered, this approximation fails at describing some models involving 
non linear laser-dipole interaction, see e.g. [17j . Consequently, the norm of the field cannot be 
always considered as a small parameter, and numerical solvers have to tolerate large controls, as 
the one described here after. 



Another distinct circumstance is when there are M systems which are exposed to the same 
laser field. Each system is characterized by its own internal Hamiltonian Hq and dipole moment 
IJLk{x)] in addition each has its own orientation denoted ^/c with respect to the incident laser 
direction. Some systems may be identical, in which case they will share the same H^ and fJ-kix) 
but may have different ^t- The governing equation is: 



(2) 



idtMx,t) = {H^{x) - ^lk{x)cos{^u)e{t))i^k{x,t), x 

0fc(x,O) = Vfco(a;), xeR^,fc = l,...,M. 



In this case the goal is to control all systems at the same time. We refer the reader to [18 1 [19 1 
120] for (positive) controllability results and numerical simulations performed with up to M = 300 
systems or even in the continuous limit ^ G [—1, 1]. 
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Throughout this paper, T > is the time of control of a quantum system. The space 
LP{0,T; X), with p G [l,+C!o) denotes the usual Lebesgue space taking its values in a Banach 
space X. The notation W^^'^(0, T) corresponds to the space of time dependent functions belong- 
ing to L^{0, T; M.) such that their first time derivative also belongs to L^{0, T). We denote by 
the space L'^[W and by W^^°° and H"^ the Sobolev spaces H/2,oo(-]^7^ r) and H'^{W, C). The 
space C{H^) is the space of linear functionals on H^. One can refer to [21^ (or the introduction 
of [22]) for more details about the definitions of these functional spaces. 

Finally, in order to introduce some numerical solver of ([T]), let us consider an integer TV and 
At > such that NAt — T. We introduce the time discretization {tj)Q<j<N of [0,T] with 
tj = jAt and we also denote by t^^i the intermediate time = (^j + ^)At. 

Let us first recall some basic results of existence and regularity of the solution of the TDSE. 
These are corollaries of a general result on time dependent Hamiltonians (see [23], p. 285, The- 
orem X.70). 

Lemma 1. Let fi e C{H^), V G W^'^ , e G L'^{0,T) and ipo G . The Schrddinger equation 

r idtm = (Ho - Mt)) m x (o, t) 

^' \^(0)=^o, M'', 

has a unique solution ip G L°°{0,T;H^) n W^'°°{0,T; L'^) such that 

IIV'(0IIl~(0,T;H2) + ||(9tV'(^)l|L~(0,T;L2) < C\\fi\\c{H'^)\\e\\w^.i(oX)Ho\\H^- 

Moreover, for all t e [0,T], ||^(t)||L2 = HVolU^. 



It is also well known (see [22] for instance) that for any T > and (f>o G H^, if we have 
e(t) = e G M, independent of time t, the Schrodinger equation 

idt^it) = (Ho ~ fie) (j){t), m X (0, T) 
m = 00, 

has a unique solution (/)(i) = S'(i)(/)o such that (f> G C([0, T]; i?^) n C^iO, T]; L^)^ where (S'(i))teM 
denotes the one-parameter semi-group generated by the operator Hq — jie. Moreover, for all 
t G [0, T], S{t) G C{H^) and we have 

S{t)^oeC{<d,T-H^), V0oGi/2; 

\\S{t)\\c(H-) <l + CT<K, me[0,T];K = KiMciH-),s^,^); 
(4) ||^(i)'/'olU^ = II^oIIl^, V^oei^VtGM; 

S'(O) = Id, 

S{t + s) = S{t)S{s), Vs,iGM. 
Therefore, the solution of Eq. ([3]) is obtained equivalently as a solution to the integral equation 

tpit) ^ S{t)'ipo + i [ S{t - s){e{s) - e)fiip{s)ds. 



3. The "toolkit" method 
Wc now present the "toolkit" method and describe the corresponding error analysis. 
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3.1. Algorithm. In this method, we assume that the control field e satisfies the following hy- 
pothesis: 

in) Vi € [0,T], e(i) e [e^i„,£^J. 

The values of the control field are discretized according to: 
(5) £f = emin + £Ae, I = 0...m, 

with m = ■ Here, the values have been here uniformly chosen in the interval 

[^inin, £max]- If somc properties of the field are known, e.g. its mean value or its variance, some 
improvement of the method can be obtained by optimizing the distribution of the values e^. 
More generally, this topic enters the field of scalar quantization, that will not be considered in 
this paper. We refer to f24| and the references therein for a review of standard methods in this 
domain. 

Remark 1. The hypothesis % often holds in practical cases. From the experimental point of 
view, there exists a technological hound for the laser amplitude. From the mathematical point of 
view, the field solves an optimality system of equations, which induces L°° bounds on the control. 
For example, consider the minimization of the functional 

J{e) = ||?A(T) - iJtargeth^ +a [ e{tfdt, 

Jo 

where a > 0, ijjtarget is a given target state and ip is the solution of ([I} (see [25] and references 
therein for details about this problem). The critical point equations read: 

( idtip{x,t) = {Ho{x) - ^i{x)e{t))ip{x,t), 
1 '>Pix,0) = ipaix), 

e{t) - --(x(i),A^V'(i))L^, 
a 

r idtx{x,t) = iHo{x) - n{x)e{t))x{x,t), 

\ X{X,T) ^i)target-'^{x,T), 

Using the above equation and thanks to the norm preservation (see one finds that 

< ^ which means that the hypothesis TL is satisfied with Smin — —■^Wf'WciH^)! 

£max = f ||M||c(_f/2)- 

In order to solve numerically equation (jS]), the "toolkit" method proceeds as follows. 

Algorithm 1. ("toolkit" method) 

(1) Preprocessing. Precompute the "toolkit", i.e. the set of propagators: 

Si{At) for£ = 0,--- ,m, 

where {Si{t))t^]gi denotes the one-parameter semi-group generated by the operator Hq — 
pL£i, the sequence {si)i=o.... ^m, being defined by ([S]). 

(2) Given a control field e e satisfuinafHl and ipQ — ipo, the sequence {ipj^)j=o,...,N that 
approximates {^{tj))j=Q^...^N, is obtained recursively by iterating the following loop: 

(a) Find: 

ij = argmin^^i ... _,„{ |e(t^+i) - e^l}, 
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(b) Setijf+,=Si^{At)i,f. 

In this "toolkit" approximation, we consider that the changes in the Hamiltonian H(t) := 
Hq — iJ,e{t) can be neglected over a time step At. In this way, if Ae = (infinite "toolkit"), and 
for a relevant time discretization, the simulation corresponding to piecewise constant control 
fields is exact (see [55] , for more details about the use of piecewise constant function in quantum 
control). Such a property does not hold with methods that approximate the exponential, e.g. 
the second order Strang operator splitting. Indeed, these approaches introduce an algebraic 
error, due to the non-commutation of the operators Hq and /z that is consequently proportional 

to ||e||L=(0,T)- 

Remark 2. In the original form of the "toolkit method" |27[ 1111 112) . the mid-point choice 
proposed in Step \2a\ of Algorithm [7] is not considered. Yet, the introduction of this strategy 
enables us to improve the order of the method (see the analysis hereafter). 

3.2. Scope and numerical considerations on the "toolkit" method. The 'toolkit' pro- 
cedure relies on the precomputation of the matrices Si{At) for multiple values of £i and is used 
in applications that require repeated resolutions of the Schrodinger equation such as control 
framework or inverse problems. The pre-computation will be a good investment as soon as the 
number of resolutions is high enough. 

However, even for the control or inverse problems not all circumstances are fitted to the use 
of the method above. The number m of matrices Si{At) to be computed is not necessarily a 
severe limitation as this can be trivially parallelized (see also our second improvement of the 
method, SecEJ which enables to greatly reduce the toolkit size). The obvious limitation arises 
when the computation of St,{At) is difficult, for instance when the system is posed in a high 
spacial dimension 7. If the matrix of Ho — /le^ in the Galerkin basis containing iV^ functions is 
not sparse the computation can scale as high as N^. 

Even then, this scaling is routine for density matrix computations: in contrast to ([l} in this 
formulation the evolving object is not the wavefunction but a density matrix operator which is 
a (self-adjoint trace-class Hilbert-Schmidt ) operator on _L^(R^)). 

It is not the object of the paper to propose novel space discretization of the multi-dimensional 
TDSE equation neither general purpose methods for computing the exponential of a matrix |28) 
so we will suppose that the user will select a setting that either avoids a full matrix Hq — /i£^ or 
manages to obtain a small Galerkin basis, such as those arising from spectral methods or reduced 
basis [29l|30] (see [31] for an application in quantum chemistry). Yet another alternative is to use 
a low rank representation of Si{At) i.e. the projection projection of Si{At) on a small number 
of eigenvectors of H^ — fiei; the eigenvectors can be computed without a full diagonalization of 
the matrix using e.g., Lanczos's method |32j . To summarize, the "toolkit" method is proposed 
in several contexts including 

- the low dimensional systems 

- any situation when a density matrix computation is possible 

- the situation in ([2]) when the problem is not an individual system by itself but a high number 
of e.g., identical systems orientated differently with respect to the laser field. Note that in this 
case the propagators Se{At) are common for various values of G [—1,1]. 

3.3. Analysis of the method. Let us now present an error analysis of the "toolkit method". 
More precisely, this section aims at proving the following result: 
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Theorem 1. Let e G W'^'°° {Q,T) and ijj the corresponding solution of ([3]). Let tp^ be the 
approximation of tp obtained with Algorithm]^ Given At > and Ae > 0, there exists Ai > 0, 
A2 > 0, that do not depend on ||e||Loo(o,T) such that: 

(6) mr) - V^^(T)|U2 < AiAe + X2At\ 
Moreover, there exists !/i > 0, 1^2 > depending on \\£\\w'-'^{o,t) such that: 

(7) IIV(T) - ^^(r)||H2 < i^iAe + v2Ae. 

Remark 3. This result shows that the "toolkit" method enables to work with large control fields, 
transferring the computational effort due to such cases to the preprocessing step: given Ae, the 
computational cost of this step only depends on the norm ||e||i,=(o,T); on Hypothesis IWl 

Proof. To obtain (|6]) and ([7]), we will focus first on the local error, i.e. the approximation 
obtained on one time step [tj,tj^i]. 

The sequence {ipf)j=o,...,N is a time discretization of the solution of 
idttp'^it) = {Ho - fie{t)) i^^it), W x (0, T) 

where the space variable has been omitted and e{t) — eg - is constant over each interval [tj,tj^i [= 
[jAt, {j + l)At[, with j = 0, - 1. We denote by {Sj{t))^^Q (instead of Sg.) the one- 

parameter semi-group generated by the operator Hq — ^ei- and we introduce 5{t) — e{t) — e 
where e (instead of e^ ) is the constant value of e{t) over [tj,tj+i\. Therefore, the solution ip of 
([3]) is actually the solution of the integral equation, settled for t G \tj,tj^i[: 

(9) %p{t) = Sj{t - tj)ijj{tj) +i [ Sj{t- s)^iS{s)ip{s) ds. 



(8) 



For the upcoming calculations, one should notice that we have 

(10) |<^(t,+i)| < ^ 
and that for all t G [tj,tj+i], 

(11) 1-^(^)1 < ^(Ae+p|U^(o,T)At). 
We consider the following decomposition: 

^(T) - V^(r) = ij{T) ~ SN-i{At)ip{tN-i) 

N-2 

+ SN-i{At) . . . S',+i(At)(^(t,+i) - S,{At)i:{t,)) 
3=0 

+ SNMAt)...So{At)^o-^''{T) 
where the last line is equal to since satisfies ^ on [0,T]. 

From now on and in all the following sections, we will consider either that UfAolli^ = 1 or that 
llV'ollff^ = 1. From @, we know that the operators Sj are isometrics in L^. Therefore, the use 
of a triangular inequality brings 

(12) MT) - ^^(T)|U. < ^ 5,(At)^(t,)IL. . 
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We will thus calculate and estimate in L'^-norm for all j the difference 





s)S{s)fiip{s) ds 






s)6{s)fi{il;{s) - 


Sj{s ~tj)iP{tj)) ds 


rtj+i 


— s)S{s)^Sj{s 


~tj)i;{tj)ds. 


Jtj 


s)S{s)fj. (V'(s) - 


Sj{s -tj)^{tj)) ds 



(13) 



S{s)(pj{s)ip{tj) ds 



where (pj{s) := Sj{tj+i — s)iJLSj{s — tj) G C{L^). 

In what follows, we work in parallel on and _ff ^-estimates of ip{tj+i) — Sj{At)ip{tj). We 
will need basic and _ff^-estimates of il;{t) — Sj{t — tj)ijj{tj) for the study of the first integral 
term of (fT3|) . the second one will be dealt with using a Taylor expansion of S{t). 

From Lemma [TJ © and (dJ), it is easy to obtain coarse estimates of ip{t) — Sj{t ~ tj)ip{tj). 
Indeed, for all t in [tj,tj^i], one can write 



\\m~Sj{t-tj)i^{t,)\\L2 



i / Sj{t ~ s)ijS{s)^jj{s) ds 



< 



tj 
tj+1 



L2 



\Sj{t - s)^lS{s)lp{s) 



1 



lL2 



ds 



< AtMciL2^- {Ae + ||e|U^(o,T)At) H^olU^ 



< 2ll/^ll£(L^) {Ae+\\eh^^o^T}At)At 



(14) 

and the i/^-cstimatc gives 

(15) U{t) - Sj{t - t,)4,{tj)\\H2 < K\\e\\w^,iio^T)\\^^\\ciH-) {Ae + ||e||i»(o^T)At) At 

where K = K{\\^\\ci^H'^),ei-na,x) is a generic constant that estimates every HS'j l|£(^f2-). 
Therefore, we can obtain more accurate estimates of the first integral term of ([T^. Thanks to 
p4p . we obtain 



Sj(tj+i ~ s)(5(s)^(-0(s) - Sj{s -tj)'4}{tj)) ds 



L2 



< 2llAi||£(L^) (A£+||£|U-(o.T)At) At sup U{s)~S,{s-t,)i;{t,)\\L2 



(16) 



< tIIa*II£(l2) (A£+||£|U^(o^T)At) At 



te[ij,tj+i] 

2 . .2 



i»(o,T)A<n At2 
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Working now on the if^-estimate, we deduce from p3)) in the same way that 

ftj+i 



(17) 



L2 



< \KM\c(m)\\e\\w^M,Q,T) (A£2 



+ l|e|li=o(o,T) 



At 



In the two cases (L^ and i?^), estimates are stronger than the ones we look for, and we can focus 
on the second integral term of we want to deal with. 



We first consider 
(18) 



s ^ Sj{tj+i - s)^iSj{s -tj) 

and note that for all "0 G H^, \\(pj{s)ip\\H2 = \\Sj{tj+i — s)iJ,Sj{s — ij)'0||^2 < iir||-0||j:^2 so that 

(19) Vs e [t,,t, + i], ||^,(s)|U(H2) < K\\fi\\ciH'^y 

Let us now consider the derivatives of <^j(s). Since {Sj{t))t<£m denotes the one-parameter semi- 
group generated by the operator Hq — /ze, the jC{H'^) identity 

holds and minor calculations give, Vs e [tj,tj+i], 

dsipjis) = iSjitj+i ~ s)[Ho,fi]Sjis -tj) 
dssf jis) = Sj{tj+i - s)[[Hq,^],Ho - fie]Sj{s - tj). 

Therefore, 

(20) \\dsV,{s)\\c^H-) < K\\[Ho,fi]\\c(H^) 
\\dss^jis)\\ciH^) < K\\[[Ho,h],Ho- fie]\\^^^^y 

If we consider the L^-analysis of the method, then (pj{s) G C{L^) and Vs £ [tj,tj+i], 

IIVj(s)ll£(L2) < ||Ai||z:(L2) 

(21) \\dsVj{s)\\ciL^^ < \\[Ho,fA\\ciL-) 
\\dssy^3is)\\c{L^) < ||[[^^0,Ai],-ffo -/i^]||£(i2) ■ 

Let us now write the third order Taylor expansion of t t-^ 6{t) = e{t) ^ e in a neighborhood 



oft 



Sis) = <5(<,+ i) + (s-i,.+ i)<5(t,.+ i) + -(s-i,.+ i)25(0(s)) 
= S{t^+i) + {.s- i )£(t,+ i) + i(s - t^+ifs{eis)), 

with 0(s) e [tj,tj-i-i]. We now focus on estimating the term i / 6{s)Lpj{s)ip{tj) ds. By 
of (pT|) and the L^-norm conservation, we obtain 



means 



S{t.^i)(pj{s)iP{tj)ds 



L2 



ANALYSIS OF THE "TOOLKIT" METHOD FOR THE TIME-DEPENDENT SCHRODINGER EQUATION 9 

rtj+i 



L2 



i t 



i'fiitj+i + s) - ^{tj+i - s))ij{tj)ds 
sdu^{u)'ipitj) duds 



t 1 — s 



L2 



^1 



< -||[i7o,/i]||£(L^)||£||L~(t„t, + i)^*' 



and 



s - < 



< ^llMll£(L2)||e||L~(t,,t,+i)A^^ 



Combining these results with (|16l) . we estimate p3)) as foUows: 

||V.(t,+i)-5,(At)V.fe)||L^ < i||M|li(i.)(A£'At' + PfooAt^) 

+ ^llMll£(L^)A£Ai 

+ ^ (2||[-ffo,A*]||£(L^)l|£||oo + ||Mll£(L^)l|e1|oo) At3, 

• lloo- By means of ([T^ . the global L^-estimate is then: 



with 



< 



T . 



T 



As 



T 
24 



;2||[-ffo,A']IU(L^)PI|oo + ||mI|£(l^)||£1|oo) At2, 



and ([S]) can be deduced with the following constants Ai and A2 independent of ||£||l°o(o,t) 
T. 



Ai 



■M\c{L-) (1 + AeA<||M||£(i2)) 



T 



T 



T 



^2 = 1 ||M|ll:(i2)||£||LAi + — ||[i/0,A']||£(L^)l|£||oc + ^MciL^Moo- 

Let us now prove the estimate. By means of (IT71) . (IT^ and (pHl) and keeping in mind that 
is a generic constant depending on ||/L(||£(i/2') and Emax, we can repeat the previous analysis 
to find the local estimate: 



+ KAeAt + K(\\[Ho,^l]\\ciH-~,\\e\\L^ + \\e\\L^)At^ 



Since one can prove that we can actually write a more precise estimate of Sj (At) and replace K 
by 1 + CAt (see properties ([1])), we get: 

||5,(Ai)l|£(H^) <l + CAt. 
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and since we have the fohowing intermediate result, where M > depends on ||/i||£(if2), Sr, 
and T but is independent of N: 

3=0 

The global estimate is obtained as follows: 

j=o 

N-1 

< 5] (1 + CAtf-^ Ikliw-(o.T) {As'Af + \\s\\lAt^) 

+ J2{l + CAtf-^ {AeAt+{\\[Ho,tA\\cim)\\i\\oo + \\e\\oo)At^[ 

< M||£||,4.i,i(o,T) {Ae^At + \\e\\i,At^) + MAe 
+ M(||[i/o,/i]||£(H2)||e||oo + ||e||oo)Ai'. 

We finally get vi and V2 and conclude the proof of Theorem [T] 
1.1 = M{1 + \\e\\w^. 1^0 ^T)T As At) 

U2 = A/(|l£||H.i,i(o,T)PllLAt+||[i/o,A^]||£(H2)||e||oo + ||e1|oo). 



□ 



Remark 4. The estimate ^ is consistent with the fact that Algorithm]^ used with a relevant 
time discretization is exact for the piecewise constant control fields. 



4. Improvement in the limit of low intensities 

We now describe a way to improve the time order of the previous algorithm. Since some 
constants in the following analysis depend in this case of the L°°-norm of the field and the 
method requires that the "toolkit" size scales At'^(einax — ^min), it applies in the case of {L°°-) 
small control fields. 

4.1. Algorithm. The algorithm we propose mixes the "toolkit" and the splitting approaches, 
in the sense that it applies sequentially various operators to correct the third order local error 
that appears in the proof of Theorem [TJ 



Algorithm 2. (Improved "toolkit" method for low intensities) 

(1) Preprocessing. Precompute the "toolkit", i.e. the set of propagators: 

Se{At) for£ = 0,--- ,m, 

where {Si{t))t<£K. denotes the one-parameter group generated by the operator Hq — fJ-S£, the 
sequence (e^)^=o,- - ,m, being defined by ([5])- Include in this set the two special elements: 
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and the initial exponents ao and (3q such that (e being extended as an even function on 
[-TM]): 

ao := = e{ti) + 0{At }, 

e{t,)-2e{t^)+e{0) 
^° ■= At^ =e{t.) + 0{At^). 

(2) Given a control field e G satisfvino[H\ and ib^^ = the sequence {^^j^)j=Q^....N 

that approximates i'4'{tj))j=o,...,N, is obtained recursively by iterating the following loop: 

(a) Find: 

= argmin^^i ... ,„{|e(ij+i/2) -ef|}, 

(b) Compute aj and jSj such that: 

(22) a, := ^^^^±^1^ = . ) + O(At^), 

e(t,+i) - 2£(t,,i) +£(<,-) 

(23) /3, := '-^^e{t^^^) + 0{At'). 

(c) Set 7/;|fi = SeJAt)n°'^Qf^^tp'^. 



In many cases, e.g. in the experimental frameworks, only the values of the field can be han- 
dled. The use of exact values for the time derivatives has then to be avoided when possible. 
This motivates the introduction of approximations (|22]) and (|23)) of e{tj) and s{tj) in the latest 
definitions. The analysis presented hereafter shows that this does not deteriorate the order of 
the method. 

In this method, one must perform two online matrices exponentiations. By working in a 
basis where one of these two matrices is diagonal, the cost of Step dc] can be reduced to one 
exponentiation, making the cost of this method equivalent the second order Strang operator 
splitting. 

4.2. Analysis of the method. We can now repeat the analysis that has been done in the proof 
of Theorem [T] to obtain the following estimate. 



Theorem 2. Let e £ W'^'°°{0,T), ip be the corresponding solution of ([3]) and ij}^^ the approxi- 
mation of ip obtained with Algorithm\^ Given At > and Ae > 0, there exists X'l > 0, A2 > 0, 
with X'l independent of ||e||L=(o,T) such that: 

\\^{T)-i;"'iT)U2 <X[Ae + \',At\ 

Proof. In the framework of this new algorithm, we note that on every time interval ]tj , tj^i [, the 
approximation ip^^ is the solution of the evolution equation: 

r tdti^'^it) =(i/o-A^e)^^^(t), W'x{t„t,+i) 
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where we set ■0(0 ) = 4'Q- We will keep the notations {Sj, S{t), (p,...) of the proof of Theorem 
[U and we first focus on the local error analysis. We consider the following decomposition: 

3=0 

+ S'Ar_i(At)f]""-ie''"-V..5o(At)f}"''e^>o-V^^^(r) 
where the last line is equal to since satisfies ([Ml) on [0,r]. 

The operators Sj are isometrics in L^, we will consider that HV'oIIl^ = 1 and we also have, for 
allj 

(25) n'^^Q^^ = g^[Ho,M]At^g^,.At^ ^ ^ (^[^0' + + ^(^*') 
and thus 

(26) Wn^'^Q^^ 11^^^,^ < 1 + C>(A<3). 
Therefore, the use of a triangular inequality brings 

(27) 1|0(T) - V'^(T)|U. < (1 + 0(At2)) 5] ||^(t,+i) - 5,(Ai)f7"^e^^-^(t,)||^, 

and we will calculate and estimate in L^-norm for all j the difference 

V'(ij+i) - Sj{M)Vt''=<df^^il;{tj) 
= 5j(Ai)V'(ij) + i / 'SjCij+i - s)5{s)^l'lp{s) ds - Sj{M)Vl°'^Q'^^'ip(tj) 



We define Y{s) = V(s) - S'j(s - tj)n°'iQ'^iil]{tj) for all s G and obtain 

(28) r(ij+i) = S'j(At) (Id - 17"^e''^) V(ii) 

+ i Sj{tj+i~ s)5{s)iiY{s)ds + i S{s)ipj{s)n°''Q^^i;{tj)ds 



where ipj{s) := Sj(ij+i — s)^Sj{s ~ tj) and its derivatives have been estimated in in (|2ip. 
As we did in Theorem [TJ we start with an estimate of the first integral term of ([^5]) . For all 
t G]tj,tj+i], we can write: 

Y{t) = m - s,{t - t^m,) + s,{t - t,)i,{t,) - s,{t - t,)rt"^eP^ij{t,) 

Sj{t - s)(5(s)^V(s) ds + Sj{t - tj) (Id - 17"^e^^) V(ii)- 
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Moreover, for all t we have 



Sj{t — s)S{s)fj,ipis) ds 



+ \\Sj{t- tj) (Id - f]"^e'^^) 



L2 



t2 



The operators Sj are isometrics in L? and Vi G [0,T], |j?/;(i)l|L2 ~ IIV'oIIl^- Therefore, we deduce 
from (US]) that 



\s,{t - 1,) (Id - i7"^eft) ^(t,)|L. < ( g ll[i?o, + ^ 



Since it is clear that wc also have 



Sj{t — s)(5(s)/iV'(s) ds 



< 



|£(L2) 



(A£+||e|U»(o,T)Af) At, 



i2 



one can finally deduce that: 
'I 

1 



Sj{tj^i — s)S{s)fj,Y{s) ds 



L2 



< -||/i||£(L2)(Ae+p|U^(o,T)At)Ai sup ||y(t)|| 



(29) < -||A*||^(i.)(Ae+||e|U..(0,T)At)'At2 + O(AeAt4) + O(Ai5). 

We focus now on the first and third terms of (|28|) . Using psp . we get 

^,(At) (Id - 17"^e^0 V^(t,) = -S,{At) (^^[Ho, m] + ^/.) V(i,)At3 + ^(t,)0(At6). 

Let us then consider the second integral term of On the one hand, we consider the fourth 

order expansion of 5 = e — e in a neighborhood of tj^i: 

6{s) = 5(t^.+ i ) + is - i )<5(t,+ i ) + lis - t^+.f6it^+.) + i(,s - t^+.f5^^\eis)) 



= Sit^^.) + is- . )£(t,+ . ) + lis - t^+^eit^+^ + \is- t^^^fe^^^Ois)) 



where 6'(s) € [tj,tj^i\. On the other hand, we calculate and/or estimate the four corresponding 
terms in 

i / (5(s)(^j(s)J7"^e'^^V(ij)<is- 
From and (HH), the term of order gives: 

1„ „o 



< 



2llA'll£(L2) 



AeAt. 
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For the term of order 1 , we can write 



At 



i i it. 



— i e [t 



iAi rt,.i+s 



t . , I —s 



sdu'fij{u)n"^e^^ip{tj) duds 
s{duVj{tj) + (u - tj)T{u))n°''e'^'-ip{tj) duds 



^2 '^3 y^Lj im,, fj,iii. Kj "vyjl 



-Sj (At) [Ha,n]n°''Q^"iJj{tj)At^ 



, 1 -s 



(m - tj)T{u - tj)n"^e^'^(tj) duds 



j^S, (At) [H^,ii\i;{t,)At' + S, {At) [H^, ti]ij[t,)0(At'^) 



lAt rt,^i+s 



s{u - tj)T{u - tj)n°'^e^'ip{tj) duds 



Jt..i~s 



where we used (f^ . (^5]) and (1^ and the function t : s G [0, At] t(s) £ C-{L^) is defined as 
the function that appears in the fohowing expansion of d^i^j around tj, for any tp £ L'^ 

du^j{u)4' = du'Pj{tj)lp + (u — tj)T(u — tj)Tp 

= iSj{At)[Ho,n]ij + {u- tj)T{u - tj)ip. 
Using the estimate (coming from ((2T|) ) 



(30) 



||t(s)||£(l2) < [Ho, fi], Ho- fie Vs G [0,At] 



along with ||V'(ij)||L2 = IIV-oIIl^ ^l,^ and ([211) we find that for all j, 



< (a,+0(Ai2)) 



lAt 



s(w - ij)r(M - tj)n°'^Qf^^'4}{tj) duds 



t - , 1 —s 
3+3 



L2 



lAt ^t,.+ l+a 



i . I 1 — s 



s(M-ij) [Ho, n], Ho - fie [[^{tj)[[^2 duds. 



- 24 



We also prove easily that for all j 



At"^ + O (At^) 



\\Sj [At) [Ho,fi]^{tj)0{Af)\\^,, = O {Af) . 

For the term of order 2, using (f23|) . ([25l) and (f2T|) and the first order expansion of ipj around tj, 

(fj{s)'4> = ipj{tj)'ip + {s- tj)eis - tj)i) for aU V' G , defining 6* : s G [0, At] i-^ e{s) G -C(L^), we 
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can write 

'3 + 1 I 



ie(t.f,i) 



24 



Sj (At) ^VfeOAt^ + Sj (At) nil;{tj)0{At^ 



{s - t^+.f{s - t,)e{s - t,W'QP^T^{t,) ds. 



Using (l2Tt . we get the estimate ||^(s)||£(i2) < Ho,fi] , Vs G [0, At], and using it with 

(j23|) and ([26|) . we obtain that for all j, 



L2 



< i(/3,+0(At2)) r^\.s-t^.+ .f(t,-s)||[Ho,M]IL(L2)||^(<,)Ili. 



< 



At 













du 



|^(^.) At^ + O {At") 



and we also prove easily that \\Sj {At) ^'\ij{tj)0{At^)\\^^ = O (At^) . 
Combining these results with ((29|) into equation (|28l) . we obtain: 

/tj+i 



< 



i2 



^^(tj+i — s)5{s)^Y {s) ds 



L2 



< i AeAt + g II [[i?o, m], - t^e 



2 

|||[i/o,/i]||£(,..)At4 + i||^|||(^. 
©(AsAi*) + ©(At^) 



At^ 



(Ae2Ai2 + ||e||i^(o.T)At^) 
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We have now a local in time estimate that should be traduced in a global one, and from (l?7l) . 
we get 

||^(r)-v^^^(r)|U. 



< 



< 



(l + 0(Ai2)) 



Af-l 

El 

J=0 



lL2 



48 



Ae- 



24 

||[i7o,M]||£(L2^At 
+ OiAeAt^) + 0(Ai4) 
The result follows, with 



[Ho,ij],Ho ~ fis 
T 



) At' + 2 IIa^II£(L^) (Ae'Ai + ||£||ioo(o,T) 



At-" 



A'i = |llA^II?:(L.)(l + AeAt) 



and 



48 



mo, 



\c{L^) + J\\^^\\c(,L^)\ 



Il°°(o,t) 



In this theorem, the constants depends on ||e||L°=(o,T) through the commutator [Hq, /i], Hq 



□ 



lie 



that appears in ((30|). This contrasts with the result obtained in Theorem[T] The explanation 
of this situation comes from the fact that the norms of tpj{s) :— Sj{tj+i—s)fiSj{s — tj) (defined in 
(fT8|) ) and its first derivative does not depend on ||£||loo(o,t), whereas its second derivative does. 
Thus, errors in Algorithm [5] depend on L°°-norm of the control field as in the case of the second 
order Strang operator splitting. Although these two methods present the same computational 
complexity, the order of Algorithm [2] is higher when Ae scales At^ . 



5. Improvement in the limit of large intensities 

We now describe a way to improve the time order of the Algorithm [T] in the case of large 
intensities. The following method enables to replace Ae by AeAt in the estimates. 

5.1. Algorithm. The algorithm we propose improve the accuracy in the approximation of e. 
This improvement is obtained by using two "toolkit" elements instead of one at each time step. 



Algorithm 3. (Improved "toolkit" method for large intensities) 

(1) Preprocessing. Precompute the "toolkit", i.e. the set of propagators: 

Si{At) for£ = 0,--- ,m, 

where {Si{t))t£m denotes the one-parameter group generated by the operator Hq — fj,e£, 
the sequence (e£)£=o,- - ,m, being defined by 

(2) Given a control field e G L°° satisfying IWl and ipo^ = V'O; the sequence {'4'j^)j=Q n 

that approximates {i^{tj))j=o,...,N, is obtained recursively by iterating the following loop: 

(a) Find £j such that: 
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(b) Compute aj and /3j such that: 

ctjSe^ + PjSi^+i = £(*j+i/2) 
(31) a,+p, = 1 

(c) Set V/fi = Si^+i{Atf^Si^{Atr^i^-J^. 

In this method, one must perform two online matrices exponentiations. The cost of the 
corresponding step, namely Step[^can be reduced to three matrix products when precomputing 
the mappings between the diagonalization basis of two consecutive "toolkit" elements. 

Remark 5. Another way to reduce the cost of this step, is to quantify the values of aj (and 
Pj) and precompute a "toolkit" containing elements of the form : S(^+i{At)^^ Se^{At)°'^ . This 
method is tested in Sec. 

5.2. Analysis of the method. We can now repeat the analysis that has been done in the proof 
of Theorem [T] to obtain the following estimate. 



Theorem 3. Let e G W^'°°{0,T), be the corresponding solution of ([3]) and ip'^^ the approxi- 
mation of ip obtained with Algorithm\^ Given At > and Ae > 0, there exists A" > 0, A2 > 0, 
both independent of ||£||l=°(o,t) such that: 

ll^(T) - i^'^'iniL^ < \'{AeAt + \'iAt\ 

Proof. In this algorithm, two control fields are involved successively in the propagation over the 
interval [tj^tj^i]. As in the previous proofs, we introduce 5{s) ~ e{s) — e(s), with 



e{s) = 

Note first that for all s e [tj^tj+i] 



£t^ s e [tj,tj + QjAi[, 
£ij+i s £ [tj + ajAt,tj+i[. 



(32) |5(s)| < A£ + ip|U^(o,T)Ai 

and denote by {Sj{t))j^g and {Sj{t)) the one-parameter semi-groups generated 

by the operators Hq — ^eg. and — fJ^Stj+i respectively. 

Following the same analysis as for Algorithm [1] we set {■ipj^)j=o n as the time discretization 
of the solution of: 



(33) 



tdt^J^it) - {Ho - fie{t))i^"'{t), X (o,r) 



where e{t) (defined right above) is constant over each interval [tj, tj -\rajAt[ and [tj -\rajAt, tj+i [, 
with j = 0, A — 1. In the same way as we obtained ([9]), the solution ip oi ^ satisfies, 

l-t j+aj At 

^{tj + ajAt) = Sj{ajAt)il>{tj) — i Sj{tj + ajAt — s)iJ.6{s)ip{s) ds 

and 

i-tj+i 

ipitj+i) = S-{PjAt)^j(tj + ajAt) - i / S'j{tj+i - s)fiS{s)i;{s) ds. 

J tj -\-aj At 
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As in (fT^ . it gives rise to: 



tj+i 



S'^{tj^i — s)fi5{s)i;{s) ds 

/tj+Qj At 

+i I Sj{l3jAt)Sj{tj + ajAt - s)nS{s)ip{s) ds 



tj+i 



S'^{tj+i - s)pl5{s) {ip{s) - S'j{s - - a.jAt)ijj{tj + ajAt)) ds 

tj+Qj At 

tj+Qj At 

S'Jl3jAt)Sj{tj + ajAt ~ s)nS{s) (V'(s) - 5^(5 - tj)i;{t.j)) ds 



(34) 



+i 



tj +aj At 

tj+QjAt 



5{s)ip'As)Sj{ajAt)ijj{tj) ds 



S'MjAt)5{s)^j{s)%l;{tj) ds. 



where <^j(s) := Sj{tj+i — s)fiSj{s ~ tj — ajAt) and ipj{s) := Sj{tj + UjAt — s)^Sj{s — tj). As 
in the proof of Theorem [1] (see right above ^2]) ) we use the appropriate decomposition 

^(T) - ^^'^{T) = V'(T) - 5Jv_i(/3jv-iAi)5jv-i(«A'-iAi)V(ijv-i) 
+ ^ 5^v_i(/3jv-iAi)^Ar-i(aiV-iAt)...5j+i(/3j+iAi)^j+i(aj+iAt) 

x(V(t,+i) - 5j(/3,At)5,(a,At)^(i,)) 

+ 5^v-i(/?A'-iA0^iv-i(aiv-iA0 . . . 5^(/3oAt)5o(«oAi)V'o - ^^"^(T) 

where the last line is equal to since ip''^ satisfies (|33| on [0,T]. We have the corresponding 
estimate (see (IT^ 'l 

JV-l 

||^(r)-^^^(T)|U. < ^ \\i;{t,+,)^Sr{P,At)S,{a,At),p{tj)\\^, 

3=0 

and we will thus calculate and estimate in L^-norm for all j the four terms of (p4)) . As in ([T6| , 
but using now the new estimate of (5, the two first terms of the right hand side of ((M)) can 
be respectively estimated by: 



(35) 



and 



(36) 



tj +Q:j At 



S'j{tj+i ~ s)n6{s) (ipis) - 5*^(5 - tj - ajAt)%l){tj + a^At)) ds 



L2 



< P,M\liL-) [As^At^ + ]^\\e\\l^^^^^^At: 



tj +CKj Ai 



S'j{f3jAt)Sj{tj + ajAt - s)fiS{s) {tP{s) - S^is - ij-)V'(ij)) 



< a," 



l£(L2) 



Ae At- + -||e|iioc(o,T) 



A<^ 
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Let us now focus on the third and fourth terms of We have: 

tj+i rtj+i 

S{s)ip'j{s)Sj{ajAt)ip{tj) ds= S{s)ip'j{tj + ajAt)Sj{ajAt)ip{tj) ds 



+ I m 



du^'j{u) du Sj{ajAt)ip{tj) ds 



'3 + 1 



5{s)S'j{PjAt)ijLSj{ajAt)%l>{tj) ds 



tj +a.j Ai 



tj+i 



5[s) I dufjiu) du Sj{ajAt)ip{tj) ds 



and 



S'0jAt)6{s)^j{s)i}{tj) ds^ S{s)S'j{(3jAt)lpj{tj + ajAt)i}{tj) ds 



tj+QjAt i-tj+oijAt 

5{s)S'^{l3jAt) i du^j{u) du ij{tj) ds 



tj -\-aj At 



6{s)S- {(3 J At) fiSj {a J At)iP{tj ) ds 



t j H-Qj At ptj +CKj Ai 

5{s)S'^{l3jAt) I du^jiu) du ip{tj) ds. 



By means of (PT|) . we have: 

5{s) ds 



tj+i r'tj+i 

e(s) - e(tj+i/2) ds+ j £{tj+i/2) - ds 
tj Jtj 



e(s) - e{tj+i/2) ds 

i 

^{0[s))-{s-tj+y2f ds, 

t-j 



where 61 (s) e [tj,tj+i]. Consequently, 
rtj+i 



(37) 



S{s)S'Jl3jAt)fiSj{ajAt)tl;{tj) ds 



L2 



< ^llA'll£(L^)l|£llL-(0,T)Ai^. 



From (I2TI) and (132]), we obtain 



(38) 



S{s) / duf'j{u) du Sj{ajAt)'tp{tj) ds 

tj+ctjAt Jtj+ajAt 



L2 



< lf3]\\[Ho,pi]\\c(L-) (^Ae + ^\\e\\L^^o,T)At^ At' 
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and similarly, we find that: 

ptj -\-aj At 
Hi 



(39) 



5{s)S'^{PjM) 



tj+a.jAt 



du^j{u) du Tpitj) ds 

1 



< ^(^jmo,^^]\\c(L^) (^A£ + -||e|U»(o,T)Ai) Ai^ 
Combining dSS]), dSS]), dSZ]), (EH) and dSH]), we obtain: 



The global estimate follows 



< 



+ ^llMll£(L^)l|e1lL~(0,T)Ai3 

+ ^ll [^^0, \\C(L-) ( + i||£||L-(0,T) At ) At'. 



N-1 

J2 - S'^{l3jAt)S,{ajAt)i:{t^)\\^, 

MlliL^) (^Ae' + l\\e\\l^^o,T)^t^^TAt 
+ ^M\cmm\L^io,T)TAf 
+^ll[^o,/i]||£(L^) ('Ae+ i||e|U^(o,T)Ai ) TAt 



and the proof of Theorem [3] is complete, with 



Ai = :^||[i?o,M]||£(L^)T+||Mll£(L2)TAe, 



A2 = ^ll[^0,A']||£(L2)l|e|U°°(0,T)T+ ^||/i||£(L2)||£|U-(0,T)T 



1 



TAt. 



□ 



6. Numerical results 

In this section, we check numerically that the order of the estimates we have obtained in this 
paper are optimal, and we compare computational costs of the methods. 

6.1. Model. In order to test the performance of the algorithms on a realistic case, a model 
already treated in the literature has been considered. The system is a molecule of HCN modeled 
as a rigid rotator. We refer the reader to [T71 [33] for numerical details concerning this system. 
As a control field, we use an arbitrary field of the form e{t) = Cmax sin(a;t), with e^a-K — 5.10"^ 
and to = 5.10"^. The parameters are chosen in accordance with usual scales considered for this 
model. The use of an analytic formula for the field enables us to work with exact values, i.e. to 
test the cases Ae = 0. 
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6.2. Orders of convergence. To test the time order, we first work with Ae = 0, with various 
values of At. The numerical orders correspond to the ones obtained in our analysis. Curves of 
convergence are depicted in Fig. [TJ 



10-' 

10-^ 

g 10"' 
E 

3 

^ 10 

I 

g 10"' 
~ 10"^ 
10"' 
10-' 

10"' 10"'' 10"' 

At/T 

Figure 1. Error with respect to At, when Ae = for "toolkit" method and 
Improved "toolkit" I method, and whenAe = cAt for Improved "toolkit" II 
method. Here, ■0""™ stands for the approximation of ip when using the "toolkit" 
method, the second order Strang operator splitting, the Improved "toolkit" I 
method and the Improved "toolkit" II. The coefficient a is the regression coef- 
ficient. 

The order with respect to Ae is also obtained numerically by using a small time step. In this 
test, the numerical order is consistent with the one obtained in Theorem [T] The convergence 
with respect to this parameter is presented in Fig. [2l 

6.3. Computational cost. In a second test, we compare the computational costs of the meth- 
ods. To do this, we look for the values of TV = and m = that enable to reach a fixed 
arbitrary error of Tol = 5.10^"^ (recall that in any case the error cannot exceed 2). For sake of 
simplicity, we only test powers of 2. In this test, we also include the quantified version of the 
Improved "toolkit" II which is described in RemarkO In our case, the parameters a and /3 were 
quantified among 100 values uniformly distributed in [0, 1]. 

These tests show that "toolkit" methods always give better results as the second order Strang 
operator splitting. 

The two improvements we propose in this paper enable to reduce respectively the global number 
of matrix products and the size of the "toolkit" , which is in agreement with the analysis we 
have done. Note that the second improvement reduce significantly preprocessing step. This fact 
makes feasible the quantified version of it, which requires intrinsically a larger "toolkit" . 



« — Splitting a=1 .9979 
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Figure 2. Error with respect to Ae, when At is smalL Here, 
the approximation of -0 when using the "toolkit" method. 



stands for 





At 


Matrix products 


m= ^ 


Strang Op. SpUtting 


16384 


32768 




Toolkit 


8192 


8192 


16384 


Improved "toolkit" I 


1024 


2048 


16384 


Improved "toolkit" II 


4096 


12288 


16 


Quantified Improved "toolkit" II 


4096 


4096 


6400 



Table 1 . Values of numerical parameters corresponding to a tolerance error of 
Tol = 5.10-3. 
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